#PGC GWAS数据分SCZ或BD去看富集的结果，原代码在 E:\代码\4-5hmc\用hg19做参考基因组——自己比对\ASM\20201216.eQTL.GWAS支持.R
setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207")
filea=read.csv("all.FDR.sig.at.least.one.add.direction.same.diff.csv",header=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.csv("at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",header=T)
file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]
allsnp=fread("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",header=T,sep="\t")
library(data.table)
allsnp=fread("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",header=T,sep="\t")
gwas.scz1=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/ckqny.scz2snpres",header=T,sep="\t")#精分的GWAS data
gwas.scz1=gwas.scz1[gwas.scz1$p<0.05,]
gwas.scz1$id=paste(gwas.scz1$hg19chrc,gwas.scz1$bp,sep=":")
gwas.scz1=data.frame(gwas.scz1$id,group.gwas=rep("scz1",dim(gwas.scz1)[1]))

gwas.scz2=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/clozuk_pgc2.meta.sumstats.ma",header=T,sep="\t")#精分的GWAS data
gwas.scz2=gwas.scz2[gwas.scz2$p<0.05,]
gwas.scz2$chr=paste("chr",gwas.scz2$chr,sep="")
gwas.scz2$id=paste(gwas.scz2$chr,gwas.scz2$pos,sep=":")
gwas.scz2=data.frame(gwas.scz2$id,group.gwas=rep("scz2",dim(gwas.scz2)[1]))

gwas.bd=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/daner_PGC_BIP32b_mds7a_0416a",header=T,sep="\t")#双相的GWAS data
gwas.bd=gwas.bd[gwas.bd$P<0.05,]
gwas.bd$CHR=paste("chr",gwas.bd$CHR,sep="")
gwas.bd$id=paste(gwas.bd$CHR,gwas.bd$BP,sep=":")
gwas.bd=data.frame(gwas.bd$id,group.gwas=rep("bd",dim(gwas.bd)[1]))

ASD=fread("./20210113/ASD.pgc.data/daner_AUT_meta14_WW_all.hg19.Mar2016_info_0.60_maf_0.05_release_Jun2017.tsv",head=F,sep="\t")
names(ASD)=c("chr","bp_hg19","snp","a1","a2","or","lb95","ub95","effect","se","p","frq_a1","nothing","info","direction")
ASD$id=paste0("chr",ASD$chr,":",ASD$bp_hg19)
ASD=ASD[ASD$p<0.05,]

ADHD=fread("./20210113/adhd_jul2017",head=T,sep="\t")
ADHD=ADHD[ADHD$P<0.05,]
ADHD$id=paste0("chr",ADHD$CHR,":",ADHD$BP)

Depr=fread("./20210113/2019.PGC.UKB.Depression.Genome-Wide.GWAS.txt",head=T,sep=" ")
Depr=Depr[Depr$P<0.05,]


names(gwas.scz1)=c("id","group.gwas.scz1")
 names(gwas.scz2)=c("id","group.gwas.scz2")
 names(gwas.bd)=c("id","group.gwas.bd")

 con=allsnp[!allsnp$id %in% as.character(filea$id),]		###remove the ASH
conid=as.character(con$id)

setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/")

gwasid =ADHD$id
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
c=length(intersect(as.character(con1$id),as.character(gwasid)))
d=dim(con1)[1]-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",c,NA,NA)

a=length(intersect(as.character(filea$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,117012-a,d),nrow = 2,byrow = T))
result[2,]=c("117K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(filea1$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,61725-a,d),nrow = 2,byrow = T))
result[3,]=c("61K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,53425-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file1$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,8544-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file2$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,807-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file3$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,200-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"./20210113/ADHD.GWAS.random.con.1.statis.csv",quote=F,row.names = F)
#	Depr
con1=read.csv("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/con_genotype.random3.anno.hg19_multianno.csv",header=T)
con1$id=paste(con1$Chr,con1$Start,sep=":")
con1=con1[!con1$id %in% filea$id,]
c=length(intersect(as.character(con1$avsnp150),as.character(Depr$MarkerName)))
d=dim(con1)[1]-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",c,NA,NA)

a=length(intersect(as.character(filea$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,117012-a,d),nrow = 2,byrow = T))
result[2,]=c("117K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(filea1$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,61725-a,d),nrow = 2,byrow = T))
result[3,]=c("61K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,53425-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file1$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,8544-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file2$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,807-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file3$avsnp150),as.character(Depr$MarkerName)))
tmp=fisher.test(matrix(c(a,c,200-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"./20210113/Depr.GWAS.random.con.3.statis.csv",quote=F,row.names = F)

###算meta.p
library(meta)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201223.GWAS")
sel=list.files(pattern = "GWAS.random.con")
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544","807","200"),
                     unique.num=c(117012,61725,53425,8544,807,200))

result=data.frame(matrix(,nrow = 6,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:6){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(rep(150000,3))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.GWAS.statis.csv",quote=F,row.names = F)